#include "headers.hh"
#include "SysUtils.h"
#include "style.h"

Int_t fgrid = -999;
Double_t fetasyscut0 = -999;
Double_t fetasyscut1 = -999;

Int_t IDx()
{
  return (fgrid/10)%5;
}

Int_t IDy()
{
  return fgrid%10;
}

Double_t GetHistX(const Int_t ib, const TH1D *hh)
{
  if(IDx()==1){
    return TMath::Sqrt(hh->GetXaxis()->GetBinLowEdge(ib)*hh->GetXaxis()->GetBinUpEdge(ib));
  }
  else if(IDx()==2){
    return hh->GetXaxis()->GetBinCenter(ib);
  }
  else{
    printf("GetHistX bad IDx %d %d\n", fgrid, IDx()); exit(1);
  }
}

//multiplication of Hist and Graph
void HistGrMul(TH1D * hraw, const TGraphAsymmErrors * gr)
{
  const Double_t x0 = gr->GetX()[0];

  Int_t rawstart = -999;
  for(Int_t ii=0; ii<hraw->GetNbinsX(); ii++){
    if(fabs(GetHistX(ii, hraw)-x0)<EPSILON){
      rawstart = ii;
      break;
    }
  }
  if(rawstart<0){
    printf("getCombFrac bad rawstart! \n");
    for(Int_t ii=0; ii<hraw->GetNbinsX(); ii++){
      printf("%d %e %e\n", ii, GetHistX(ii,hraw), x0);
    }
    exit(1);
  }

  for(Int_t ip=0; ip<gr->GetN(); ip++){
    const Double_t xx = gr->GetX()[ip];
    
    const Int_t iraw = rawstart+ip;
    
    if(fabs(xx- GetHistX(iraw, hraw))>EPSILON){
      printf("wrong gr x ip %d x0 %e %e", ip, x0, GetHistX(iraw, hraw)); exit(1);
    }
    
    const Double_t yy = hraw->GetBinContent(iraw);
    const Double_t gg = gr->GetY()[ip];
    hraw->SetBinContent(iraw, yy*gg);
  }
}

//for fraction, i.e. IDy()==0, hraw is total yield, gr is fraction
//for to-pion ratio, i.e. IDy()==3, hraw is pion yield, gr is to-pion ratio
void getCombFrac(const Int_t nvar, TH1D * hraw[], const TGraphAsymmErrors * gr[], const Bool_t kstat)
{
  //Int_t GetCommonLimits(const Int_t ngr, const TGraphAsymmErrors** grs, Int_t* id0s, Int_t* id1s)
  Int_t id0s[200];
  Int_t id1s[200];
  const Int_t np = SysUtils::GetCommonLimits(nvar, gr, id0s, id1s);
  if(np==0){
    printf("wrong common limits !\n"); exit(1);
  }

  //for both z and pT
  const Double_t x0 = gr[0]->GetX()[id0s[0]];

  //========================================================================
  //align
  Int_t rawstart = -999;
  for(Int_t ii=0; ii<hraw[0]->GetNbinsX(); ii++){
    if(fabs(GetHistX(ii, hraw[0])-x0)<EPSILON){
      rawstart = ii;
      break;
    }
  }
  if(rawstart<0){
    printf("getCombFrac bad rawstart! \n"); 
    for(Int_t ii=0; ii<hraw[0]->GetNbinsX(); ii++){
      printf("%d %e %e\n", ii, GetHistX(ii,hraw[0]), x0);
    }
    exit(1);
  }

  //========================================================================
  TGraphAsymmErrors * gout= new TGraphAsymmErrors;
  gout->SetName(gr[0]->GetName());

  gr[nvar] = gout;

  for(Int_t icom=0; icom<np; icom++){

    const Int_t i0 = id0s[0]+icom;   
    const Double_t xx = gr[0]->GetX()[i0];

    const Int_t iraw = rawstart+icom;

    Double_t totyield = 0;
    Double_t paryield = 0;

    Double_t sumw2l = 0;
    Double_t sumw2h = 0;

    Double_t sumw1l = 0;
    Double_t sumw1h = 0;
    
    Double_t ymin = 1e10;
    Double_t ymax = -1e10;

    for(Int_t ivar=0; ivar<nvar; ivar++){
      const Int_t ip = id0s[ivar]+icom;
    
      if(fabs(xx- gr[ivar]->GetX()[ip])>EPSILON){
        printf("wrong gr x ip %d ivar %d x0 %e %e", ip, ivar, x0, gr[ivar]->GetX()[ip]); exit(1);
      }
      if(fabs(xx- GetHistX(iraw, hraw[ivar]))>EPSILON){
        printf("wrong gr x ip %d ivar %d x0 %e %e", ip, ivar, x0, GetHistX(iraw, hraw[ivar])); exit(1);
      }
    
      //totyield, or pion yield
      totyield += hraw[ivar]->GetBinContent(iraw);
      paryield += hraw[ivar]->GetBinContent(iraw)*gr[ivar]->GetY()[ip];

      sumw2l += TMath::Power(hraw[ivar]->GetBinContent(iraw)*gr[ivar]->GetEYlow()[ip] ,2);
      sumw2h += TMath::Power(hraw[ivar]->GetBinContent(iraw)*gr[ivar]->GetEYhigh()[ip] ,2);

      sumw1l += TMath::Abs(gr[ivar]->GetEYlow()[ip]);
      sumw1h += TMath::Abs(gr[ivar]->GetEYhigh()[ip]);

      if(gr[ivar]->GetY()[ip]>ymax){
        ymax = gr[ivar]->GetY()[ip];
      }
    
      if(gr[ivar]->GetY()[ip]<ymin){
        ymin = gr[ivar]->GetY()[ip];
      }
    }

    const Double_t frac = paryield/(totyield+EPSILON);
    Double_t el = -999, eh= -999;

    Double_t globalvar = 0;
    for(Int_t ivar=0; ivar<nvar; ivar++){
      const Int_t ip = id0s[ivar]+icom;

      //difference between mean and ivar
      const Double_t localvar = gr[ivar]->GetY()[ip] - frac;
      //weight to apply on deltay
      const Double_t localweight = hraw[ivar]->GetBinContent(iraw)/(totyield+EPSILON);

      //linear
      //globalvar += localweight * TMath::Abs(localvar);

      //quadratic
      globalvar += TMath::Power(localweight * localvar, 2);
    }

    //for quadratic sys comb
    globalvar = TMath::Sqrt(globalvar);

    if(kstat){
      //this overwrites all errors
      /*
      if(IDy()==0){
        el = eh = TMath::Sqrt(frac*(1-frac)/(totyield+EPSILON));
      }
      else if(IDy()==3){
        el = eh = frac*TMath::Sqrt(1/(totyield+EPSILON) + 1/(paryield+EPSILON)); //correct since in this internal conversion totyield and paryield are both real counts
      }
      else{
        printf("wrong IDy %d\n", IDy()); exit(1);
      }
      */

      //checked -> good, perfect PID = MC truth -> use it!
      el = TMath::Sqrt(sumw2l)/(totyield+EPSILON);
      eh = TMath::Sqrt(sumw2h)/(totyield+EPSILON);
    }
    else{
      //this is the long existing recipe -- simple mean of sys. error, big enough good enough
      el = sumw1l/nvar;
      eh = sumw1h/nvar;

      //add eta variation for all pion, electron and high pt proton kaon
      //it has to be in!
      if(xx > fetasyscut0 && xx < fetasyscut1 ){
        //const Double_t sysh = ymax - frac;
        //const Double_t sysl = frac - ymin;

        //maybe too large testing alternative below
        //Maxmin
        //const Double_t sysh = 0.5*(ymax-ymin);
        //const Double_t sysl = 0.5*(ymax-ymin);

        //Quadratic or linear
        const Double_t sysh = globalvar;
        const Double_t sysl = globalvar;
        
        el = TMath::Sqrt(el*el+sysl*sysl);
        eh = TMath::Sqrt(eh*eh+sysh*sysh);
        printf("testtest additional xx %e sysh %e sysl %e\n", xx, sysh, sysl);
      }
    }

    gout->SetPoint(icom, xx, frac);
    gout->SetPointError(icom, gr[0]->GetEXlow()[i0], gr[0]->GetEXhigh()[i0], el, eh);
  }
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=1){
    printf("argc!=1\n"); return 1;
  }

  const TString pwd=gSystem->pwd();

  TString vars[100];
  Int_t nvar=-999;

  Double_t jpt0 = -999;
  Bool_t keta = kFALSE;

  if(pwd.Contains("recom510")){
    vars[0]="JPT5-6";
    vars[1]="JPT6-7";
    vars[2]="JPT7-8";
    vars[3]="JPT8-9";
    vars[4]="JPT9-10";
    vars[5]="JPT5-10";//only for output
    nvar = 5;

    keta = kFALSE;
  }
  else if(pwd.Contains("recomEta")){
    if(pwd.Contains("JPT5-6")){
      jpt0 = 5;
    }
    else if(pwd.Contains("JPT5-7")){
      jpt0 = 5;
    }
    else if(pwd.Contains("JPT6-7")){
      jpt0 = 6;
    }
    else if(pwd.Contains("JPT7-8")){
      jpt0 = 7;
    }
    else if(pwd.Contains("JPT8-9")){
      jpt0 = 8;
    }
    else if(pwd.Contains("JPT9-10")){
      jpt0 = 9;
    }
    else if(pwd.Contains("JPT5-10")){
      jpt0 = 5;
    }
    else if(pwd.Contains("JPT10-15")){
      jpt0 = 10;
    }
    else if(pwd.Contains("JPT10-16")){
      jpt0 = 10;
    }
    else if(pwd.Contains("JPT10-17")){
      jpt0 = 10;
    }
    else if(pwd.Contains("JPT10-20")){
      jpt0 = 10;
    }
    else if(pwd.Contains("JPT15-20")){
      jpt0 = 15;
    }
    else{
      printf("bad pwd for jpt0 %s\n", pwd.Data()); 
      //exit(1);
      printf("should only be critical for z, i.e. IDx()==2\n");
    }

    nvar=-999;
    //pPb
    if(!pwd.Contains("JPT")){
      if(pwd.Contains("IDFF")){
        vars[0]="ETA0";
        vars[1]="ETA1";
        vars[2]="ETA2";
        vars[3]="EUA0";
        vars[4]="EUA1";
        vars[5]="EUA2";
        vars[6]="ETAall";//only for output
        nvar=6;
      }
      else if(pwd.Contains("Incl")){
        vars[0]="ETA0";
        vars[1]="ETA1";
        vars[2]="ETA2";
        vars[3]="ETA3";
        vars[4]="EUA0";
        vars[5]="EUA1";
        vars[6]="EUA2";
        vars[7]="EUA3";
        vars[8]="ETAall";
        nvar=8;
      }
    }
    //pp //for JPT (i.e. PP) there is no EUA
    else{ 
      vars[0]="ETA0";
      vars[1]="ETA1";
      vars[2]="ETA2";
      vars[3]="ETAall";//only for output
      nvar=3;
    }

    if(nvar<0){
      printf("getCombFrac nvar<0 %s %d\n", pwd.Data(), nvar); exit(1);
    }

    keta = kTRUE;

  }
  else{
    printf("bad pwd for var %s\n", pwd.Data()); exit(1);
  }

  if(pwd.Contains("GrID10")){
    fgrid = 10;
  }
  else if(pwd.Contains("GrID20")){
    fgrid = 20;
  }
  else if(pwd.Contains("GrID13")){
    fgrid = 13;
  }
  else if(pwd.Contains("GrID23")){
    fgrid = 23;
  }
  else if(pwd.Contains("GrID60")){
    fgrid = 60;
  }
  else if(pwd.Contains("GrID70")){
    fgrid = 70;
  }
  else if(pwd.Contains("GrID63")){
    fgrid = 63;
  }
  else if(pwd.Contains("GrID73")){
    fgrid = 73;
  }
  else{
    printf("bad pwd for id %s\n", pwd.Data()); exit(1);
  }

  //==========================================  

  const Int_t ntype = 4;
 
  TList * ll = new TList;

  TH1D * pionraw[nvar];
  for(Int_t ivar = 0; ivar<nvar; ivar++){
    pionraw[ivar] = 0x0;
  }

  const Int_t parid[]={1,0,2,3};
  for(Int_t pid=0; pid<ntype; pid++){
    const Int_t itype = parid[pid];

    const TGraphAsymmErrors * fracstat[nvar+1];
    const TGraphAsymmErrors * fracsys[nvar+1];

    const TGraphAsymmErrors * ratiostat[nvar+1];
    const TGraphAsymmErrors * ratiosys[nvar+1];
    
    TH1D * hraw[nvar+1];
    for(Int_t itmpv = 0; itmpv<=nvar; itmpv++){
      fracstat[itmpv] = fracsys[itmpv] = ratiostat[itmpv] = ratiosys[itmpv] = 0x0;
      hraw[itmpv] = 0x0;
    }

    Int_t nread = 0;
    for(Int_t itmpv = 0; itmpv<nvar; itmpv++){

      TFile * fraw = new TFile(Form("raw%s.root", vars[itmpv].Data()));
      if(!fraw->IsOpen()){
        printf("no raw file! %s\n", fraw->GetName()); exit(1);
      }

      TFile * fracfile = new TFile(Form("frac%s.root", vars[itmpv].Data()));
      TFile * ratiofile = 0x0;

      if(!fracfile->IsOpen()){
        printf("no fraction file %s\n", fracfile->GetName()); 
        exit(1);
        //continue;
      }
      if(IDy()==3){
        ratiofile = new TFile(Form("ratio%s.root", vars[itmpv].Data()));
        if(!ratiofile->IsOpen()){
          printf("no ratiotion file %s\n", ratiofile->GetName()); 
          exit(1);
          //continue;
        }
      }

      //====================

      //remember to increase nread in the end of loop
      fracstat[nread] = 0x0;
      fracsys[nread] = 0x0;

      ratiostat[nread] = 0x0;
      ratiosys[nread] = 0x0;

      hraw[nread] = 0x0;
        
      fracstat[nread]=(TGraphAsymmErrors*)fracfile->Get(Form("grstat%d",itype));
      if(!fracstat[nread]){
        printf("no fracstat nread %d itype %d in %s\n", nread, itype, fracfile->GetName()); fracfile->ls(); exit(1);
      }
      ll->Add(fracstat[nread]->Clone(Form("%sfile%dpar%d",fracstat[nread]->GetName(), nread, itype)));
      
      fracsys[nread]=(TGraphAsymmErrors*)fracfile->Get(Form("grsys%d",itype));
      if(!fracsys[nread]){
        printf("no fracsys nread %d itype %d in %s\n", nread, itype, fracfile->GetName()); fracfile->ls(); exit(1);
      }
      ll->Add(fracsys[nread]->Clone(Form("%sfile%dpar%d",fracsys[nread]->GetName(), nread, itype)));
        
      //==========

      TH2D * hvar=0x0;
      if(IDx()==1){
        hvar=(TH2D*) fraw->Get("kk");
      }
      else if(IDx()==2){
        hvar=(TH2D*) fraw->Get("zz");
      }
      else{
        printf("bad IDx %d %d\n", IDx(), fgrid); exit(1);
      }
      
      if(!hvar){
        printf("no hvar in %s\n", fraw->GetName()); exit(1);
      }
      ll->Add(hvar->Clone(Form("%sfile%dpar%d",hvar->GetName(), nread, itype)));

      hraw[nread] = hvar->ProjectionY(Form("raw%d", nread));
      
      if(IDx()==1) {
        style::ToNaturalScale(hraw[nread]);
      }

      //================    
      if(IDy()==3){
        ratiostat[nread]=(TGraphAsymmErrors*)ratiofile->Get(Form("grstat%d",itype));
        if(!ratiostat[nread]){
          printf("no ratiostat nread %d itype %d in %s\n", nread, itype, ratiofile->GetName()); ratiofile->ls(); exit(1);
        }
        ll->Add(ratiostat[nread]->Clone(Form("%sfile%dpar%d",ratiostat[nread]->GetName(), nread, itype)));

        ratiosys[nread]=(TGraphAsymmErrors*)ratiofile->Get(Form("grsys%d",itype));
        if(!ratiosys[nread]){
          printf("no ratiosys nread %d itype %d in %s\n", nread, itype, ratiofile->GetName()); ratiofile->ls(); exit(1);
        }
        ll->Add(ratiosys[nread]->Clone(Form("%sfile%dpar%d",ratiosys[nread]->GetName(), nread, itype)));
      }

      //=====================
      if(itype==1 && IDy()==3){
        pionraw[nread] = (TH1D*) hraw[nread]->Clone(Form("pionraw%d",nread));
        HistGrMul(pionraw[nread], fracstat[nread]);
      }

      //good!!
      nread++;
    }

    if(keta){
      //default apply everywhere
      fetasyscut0 = -1e10;
      fetasyscut1 =  1e10;
      
      if(itype==0){
        fetasyscut0 = 1;
        fetasyscut1 = 1e10;
      }
      else if(itype==2){
        fetasyscut0 = 0.5;
        fetasyscut1 = 1e10;
      }
      
      if(IDx()==2){
        fetasyscut0/=jpt0;
        fetasyscut1/=jpt0;
      }
    }
    else{
      //default apply no where
      fetasyscut0 =  1e10;
      fetasyscut1 = -1e10;
    }

    printf("testtest itype %d fetasyscut %e %e\n", itype, fetasyscut0, fetasyscut1);

    TH1D ** hin                  = IDy()==3? pionraw   : hraw;
    const TGraphAsymmErrors ** gstatin = IDy()==3? ratiostat : fracstat;
    const TGraphAsymmErrors ** gsysin  = IDy()==3? ratiosys  : fracsys;

    gstatin[nread]=0x0;
    getCombFrac(nread, hin, gstatin, kTRUE);
    gsysin[nread]=0x0;
    getCombFrac(nread, hin, gsysin,  kFALSE);

    ll->Add((TGraphAsymmErrors*)gstatin[nread]);
    ll->Add((TGraphAsymmErrors*)gsysin[nread]);
  }
  
  TFile *fout = new TFile(Form("Combined_%s.root", vars[nvar].Data()),"recreate");
  ll->Write();
  ll->ls();
  fout->Save();
  fout->Close();

  printf("getCombFracdonedone!\n");

  return 0;
}
